General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 





> 



LU 








DEPARTMENT OP MECHANICAL ENGINEERING AND MECHANICS 
SCHOOL OP ENGINEERI>:C 
OLD DOMINION UNIVERSITY 
NORPOLR, VIRGINIA 


APPLICATION OP THE METHOD OP LINES FOR SOLUTIONS 
OP THE NAV7.ER-STOITES EQUATIONS USING J NONUNIFORM 
GRID DISTRIBUTION 


By 

Jamchid S. Abolha>«ani 
and 

S. N. Tiwari, Principal Investigator 


Progress Report 

For the period ending June 30, 1983 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 



Under 

Research Grant NCCI-*68 

Robert E. Smith, Jr., Technical Monitor 


If 

V 


A 

< 



(NAS A-Ca- 173180) APPLICATION OF THE flETHOD 
OF LINES FOB SOLUTIONS OF THE N AV I EB-STOK ES 
EQUATIONS USING A NONUNIFOBH GHID 
CISTfilBUTION Progress Report, period ending 
30 Jun. '983 (Old Dominion U 


G3/6 4 


N84-I6870 ; 


Unclas 

15044 


August 1983 




DEPARmNT OF MECHANICAL ENGINEERING AND MECHANICS 
SCHOOL OF ENGINEERING 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 


APPLICATION OF THE METHOD OF LINES FOR SOLUTIONS 
OF THE NAVIER-STOKES EQUATIONS USING A NONUNIFORM 
GRID DISTRIBUTION 


By 

Jamshid S. Abolhassani 
and 

S. N. Tiwari, Principal Investigator 


Progress Report 

For the period ending June 30, 1983 


Prepared for the 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Research Grant NCCI-68 

Robert E. Smith, Jr.* Technical Monitor 


4 ^ 


Submitted by the 

Old Dominion University Research Foundation 

P.O. Box 6369 

Norfolk, Virginia 23508 


August 1983 


FORWARD 


This is a progress report on the work completed on the reseR;.-:'? project 
"I^merical Solution of Axisymmetric Flow with Arbitrary Wall Geometry." The 
period of performance on this research was June 1, 1982 through June 30, 
1983. The work was supported by the NASA/Langley Research Center through 
the Cooperative Agreement NCCl-68. The cooperative agreement was monitored 
by Dr. Robert E. Smith, Jr., of the Analysis and Computational Division 
(Computer Science, and Applications Branch), MS/12S. 


ii 


TABLE OF CONTENTS 


Page 

FORWARD ii 

LIST OF FIGURES v 

LIST OF SYMBOLS viii 

SUMMARY X 

1. INTRODUCTION 1 

2. THEORETICAL FORMULATION 5 

2.1 Governing Equations.. 3 

2.2 Stream Function-Vorticity Formulation 6 

2.3 Coordinate Transformation 8 

3. COMPUTATIONAL PROCEDURES 15 


3.1 Stream Function Equation 15 ^ 

3.2 Application of Successive Over Relaxation 17 ] 

3.3 The Application of the Method of Lines i 

to the Vbrticity Equation 20 < 

3.3.1 Advantages and Disadvantages of the Method 

of Lines 21 \ 

I 

3.3.2 Convergence, Accuracy, and Stability..... 21 | 


4. STIFFNESS ANALYSIS 

4.1 Effect of Grid Concentration on Stiffness.. 

4.2 Effect of Differencing Schemes on Stiffness 

5. BOUNDARY AND INITIAL CONDITIONS 

5.1 Inlet Condition 

5.2 Outlet Condition 

5.3 Wall Condition 

5.4 Symmetry Condition 

5.5 Moving Wall Condition... 

5.6 Initial Condition 

6. GRID GENERATION TECHNIQUES 

6.1 Algebraic Method 

6.2 Differential Method.,... 

7. PHYSICAL APPLICATIONS 

7 . 1 Internal Flows 

7.2 Driven Cavity... 


Ill 


23 

24 
27 

36 

37 
39 

41 

42 
42 

44 

45 

46 
48 

50 

51 
53 


* - 



TABLE OF CONTENTS - CONCLUDED 


Page 

8. RESULTS AND DISCUSSION 56 

9. CONCLUSIONS 92 

REFERENCES.... 95 


IV 


LIST OF FIGURES 


F igure Page 

2.1 Coordinate systems.. 9 

3.1 Node numbering scheme 18 

4.1 Stiffness characteristic for heat equation 28 

4.2 Stiffness characteristic for fluid flow equation........ 31 

4.3 Stability region for Adams-Bashforth (first order) 32 

4.4 Variation of step size with Reynolds nunber and 

number of grids point 34 

5.1 Boundary conditions for various geometries 38 

7.1 Physical model for duct flows 32 

7.2 Curved-wall diffuser geometry 34 

8.1 Grid distribution for parallel plates and pipe... 37 

8.2a Velocity distribution for flow between parallel 

plates. Re ~ 200 38 

8.2b Stream function contours for flow between parallel 

plates, Re = 200 39 

8.2c Vorticity contours for flow between parallel 

plates. Re =* 200 60 

8.3a Velocity distribution for pipe flow. Re = 200 61 

8.3b Stream function contours for pipe flow. Re = 200 6 2 

8.3c Vorticity contours for pipe flow. Re = 200 63 

8.4a Velocity distribution for flow between parallel 

plates. Re = 1000 64 

8.4b Stream function contours for flow between parallel 

plates. Re - 1000 63 

8.4c Vorticity contours for flow between parallel 

plates. Re - 1000 66 

8.5a Velocity distribution for pipe flow. Re = 1000 67 

8.5b Stream function contours for pipe flow. Re =* 1000 68 





II *. 

■A* . 


if- 




ft'.. 

i ^ t* 


|Bi 


■ I 


V 


LIST OF FIGURES - CONTINUED 

Figure Page 

8.5c Vorticity contours for pipe flow, Re “ 1000 69 

8.6 Grid distribution for a diffuser, 5=1, 9=5, 

L » 20, Lj. = 32 71 

8.7a Velocity distribution for a bell-type diffuser, 

a - 1, 0-5, L - 20, Lt - 32, Re » 200 72 

8.7b Vorticity contours for a bell-type diffuser, 

a » 1, e « 5, L=20, Lt “32, Re = 200 73 

8.7c ^tream function contours for a bell-type diffuser, 

a » 1, 9=5, L - 20, L^. = 32, Re = 200 74 

8.8 £rid distribution for a bell-type diffuser, 

a » 1, 0-10, L = 20, » 40 75 

8.9a Velocity distribution for a bell-type diffuser, 

a » 1, 9-10, L = 20, =40, Re = 200 76 

8.9b ^trean function contours for a bell-type diffuser, 

a » 1, 0 = 10, L = 20, Lg -40, Re = 200 77 

8,9c \^orticity contours for a bell-type diffuser, 

a » 1, 9 = 10, L = 20, * 40, Re = 200 78 

8.10 Uniform grid distribution for a driven cavity... 79 

8.11a Stream function distribution for uniform grid 

distribution. Re = 100 80 

8.11b Vorticity distribution for uniform grid 

distribution. Re = 100 81 

8.12a Stream function distribution for uniform grid 

distribution. Re = 1000 82 

8.12b Vorticity distribution for uniform grid dis- 
tribution, Re = 1000... 83 

8.13a Stream function distribution for uniform grid 

distribution. Re = 10,000 84 

v± 


LIST OF FIGURES - CONCLUDED 


Figure Page 


8»13b VorCicity distribution for uniform grid distri- 
bution, Re " 10,000 85 

8.14 Nonuniform grid distribution for a driven cavity....... 87 

8.15a Streau function distribution for nonuniform grids. 

Re » 100 88 

8.15b Vorticity distribution for nonunifoim grids. Re “ 100.. 89 

8.16a Strean function distribution for nonuniform grids. 

Re = 1,000 90 

8.16b Vorticity distribution for nonuniform grids. 

Re - 1,000 91 


VI 1 


LIST OF Slft.lBOLS 


a, b, c 

AI, AO 

Al, 31 
A2, B2 

H(X) 

IM 

|J| 

k 

K 

L 

NI, NO 


P 

r 

^2 

Re 

RE 

T 

t: 



element of I-iatrix (A), Eq. (4.13) 

coefficient of Eqs. .1) and (5.5) 

inflow and outflow parameter, Eqs. (5.1) and 
(5.5) 

coefficient of Eq. (2.21) 

coefficient of Eq. (2.22) 

blending function, Eqs. (6.11) and (6.12) 

step function, Eq. (7.1) 

imaginary part of the eigenvalue 

Jacobian of transformation 

iteration number 

concentration factor, Eq. (2.3) 

characteristic length 

total length of diffuser 

inflow and outflow parameter 

normalized pressure 

radial coordinate (physical) 

diffuser outlet radius 

Reynolds number 

real part of the eigenvalue 

temperature 

non-dimensional time 

non-dimensional velocity components 

free stream velocity 


viii 


LIST OF SYMBOLS - CONCLUDED 


U 


max 


z 

a 

a 


n 

6 

0 

P 

U 

M 00 
'P 

Hi 

“o 

X 


Subscripts 


maximun outlet velocity 
axial coordinate (physical) 
wall contouring parameter 
M/Re 

radial (vertical) coordinate (computational) 

index, zero for two dimensional, one for 
axis ymme trie flow 

angle of diffuser 

density 

non-dimensional viscosity 

free stream viscosity 

strean function 

vorticity 

relaxation factor 

axial coordinate (computational) 

eigenvalue 


i. j. N 


node number 


APPLICATION OF THE METHOD OF LINES FOR SOLUTIONS OF 
THE NAVIER-STOKES EQUATIONS USING A NONL7UFORM GRID DISTRIBUTION 

By 

J. S. Abolhassani^ and S» N. Tiwari^ 

SUMMARY 

The feasibility of the method of lines is investigated for solu- 
tions of physical problems requiring nonuniform grid distributions. To 
attain this, it was also necessary to investigate the stiffness charac- 
teristics of the pertinent equations. For specific applications, the 
governing equations considered are those for viscous. Incompressible, 
two-dimensional and axisymmetric flows. These equations are transform- 
ed from the physical domain having a variable mesh to a computational, 
domain with a uniform mesh. The two governing partial differential 
equations are the vorticity and stream function equations. The method 
of lines is used to solve the vorticity equation and the successive 
over relaxation technique is used to solve the stream function equa- 
tion. 

The method is applied to three laminar flow problems. These are; 
the flow in ducts, curved-wall diffusers, and a driven cavity. Results 
obtained for different flow conditions are in good agreement with 
available analytical and numerical solutions. The viability and 
validity of the method of lines are demonstrated by its application to 
Navier-Stokes equations in the physical domain having a variable mesh. 

^ Graduate Research Assistant, Department of Mechanical Engineering and, 
Mechanics, Old Dominion University, Norfolk, Virginia 23508. 

Eminent Professor, Department of Mechanical Engineering and Mechanics 
Old Dominion University, Norfolk, Virginia 23508. 


Chapter 1 
INTRODUCTION 

In engineerii'i^ and sciences, most physical phenomena may be de- 
scribed by a set of differential equations and boundary condition equa- 
tions. These equations are mostly nonlinear in nature, and in a major- 
ity of cases, they can be solved only by numerical approaches. In 
particular, the field of fluid dynamics is governed by a specific set 
of nonlinear partial differential equations called the Navier-Stokes 
equations. In certain cases, flow fields may be described accurately 
by the incompressible form of the Navier-Stokes equations* For the 
present study, the incompressible form of the Navier-Stokes equations 
is expressed in terms of the stream function and vorticity. The 
resulting equations are two coupled nonlinear partial differential 
equations. The boundary conditions for the dependent variables, stream 
function and vorticity, are expressed by a set of coupled linear 
equations which depend on the nature of physical applications. 

In the field of computational fluid dynamics, finite difference 
and finite element methods have a well developed history. Another 
method called "The Method of Lines (MOL)" has also received special 
attention for the numerical solution of certain partial differential 
equations. Detailed discussions on the method of lines are available 
in [1-8]*. If the governing equations contain both time and space 
variables, the procedure is to discretize the space variable components 
and treat the time varying component as a continuum. This leads to a 
system of coupled ordinary differential equations (ODE'S) which can be 


*The number in brackets indicate references . 


integrated with sophisticated ordinary differential equation software 
[9] having automated time-step control. The advantages of this 
approach are that the techniques can be applied quickly and the magni- 
tude of knowledge about ordinary differential equation solvers can be 
utilized for obtaining specific solutions, 

The early development of the method of lines was in the Soviet 
Union. Liskovet’s article [1] is a review of 154 papers which date up 
to the mid-sixties. This review has considered linear partial differ- 
ential equations of elliptici parabolic, and hyperbolic types. Laser 
and Harrison [2], and Hicks and Wei [3] showed extensively the 
viability and validity of the method for linear partial differential 
equations. Later, Klunker et al. [4] used the method of lines to 
calculate nonlinear conical flows. It was concluded that the method of 
lines is a useful and versatile procedure for structuring the numerical 
^.'/iutions to nonlinear equations. Jones et al. [5] presented an 
extensive discussion for application of the method of lines in elliptic 
systems. It was concluded that a large nuiiber of lines may cause 
nonconvergence. However, this conclusion was based on results obtained 
from linear systems. Loeb and Schiesaer [10] presented an elegant way 
to analyze the stability of the method of lines. It was concluded that 
higher order finite difference approximations to the spatial deriva- 
tives would improve the accuracy, stability and computational cost. 
Madsen and Sincovec [6] applied the method of lines for solution of 
several nonlinear partial differential equations. These problems were 

the diffusion of electrolytes, flow through porous media, and global 
atmospheric transport with kinetics. It was concluded that the method 
of lines gives satisfactory and reliable results which could not be 


2 


obtained by using the finite difference schemes. Hamilton [7] obtained 
solutions of axisymmetric and two'-dimensional inviscid flow over blunt 
bodies using the method of lines and observed very accurate solutions 
using just a few lines. Kurtz et al. [8] applied the method of lines 
to the viscous strean function - vorticity equations in a rectangular 
coordinate system, the particular problem discussed was the flow in a 
driven cavity. It was concluded that the method of lines is capable of 
producing solutions to the strean function - vorticity equation at very 
high Reynolds numbers where standard finite difference techniques fail 
[ 11 ]. 

The literature survey indicates that the method of lines has been 
applied by several investigators in numerical experimentations using 
only the uniform grid distribution. For many physical problems, 
however, it becomes essential to have irregular grid distributions . 

Tnis may be due to the complexity of the physical boundary geometry, 
and/or local grid resolution. Therefore, there exists a strong need 
for investigating the feasibility of the method of lines for physical 
domains that are covered with variable grids which conform to the 
boundary contours and may be concentrated in specified regions. 

The objective of this study is to establish the validity and via- 
bility of the method of lines where there is an arbitrary grid distri- 
bution. Also, this study investigates the effects of grid concentra- 
tion, Reynolds number, boundary conditions, and differencing schemes on 
stiffness and stability by using a one-dimens ions ^ advection - diffu- 
sion fluid flow equation and heat equation. 
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For viscous incompressible flow, the equations of raoti!r.i are 
derived in two-dimensional and axis >mme trie coordinate systems in 
Chap. 2. These equations are transformed from the physical domain 
having a variable grid to a computational domain with a uniform grid. 
The resulting equations are solved numerically. The method of lines is 
used to solve the vorticity equation and successive overrelaxation is 
used to solve the stream function equation. The computation procedure 
is presented in Chap. 3. The stiffness analysis is presented in 
Chap. f\: A discussion on appropriate boundary and initial conditions 

is given in Chap. 5. Information on grid generation is presented in 
Chap. 6. In this study, the grids are generated by an algebraic method 
which transforms the irregular physical domain into a uniform 
computational domain. For physical applications, specific problems 
considered are: the flow between horizontal ducts, curved-wall diffus- 

ers, and flow in a driven cavity. The applications are described in 
Chap. 7. Results are obtained for different flow conditions and are 
compared with available analytical and mmerical solutions in Chap. 8. 
The viability and validity of the method of lines are illustrated by 
its applications to the incompressible Navier-S'eokes equations. 
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Chapter 2 

THEORETICAL FORMULATION 

The theory of fluid dynamics is based upon a set of governing 
equations called the Navier-Stokes equations. For multi-dimensional 
flow* the equations are second order nonlinear parabolic-elliptic 
partial differential equations. The first-order boundary-layer forms 
of the Navier-Stokes equations are parabolic in nature and offer some 
computational conveniences. Thes^ equations* however* are not 
applicable in many realistic flow conditions such as reverse and 
separated flows. Therefore* it becomes essential to make use of the 
full Navier-Stokes equations in many flow situations of practical 
interest. A brief discussion of the basic governing equations used in 
this study is presented in this chapter. 

2.1 Governing Equations 

For viscous incompressible flow the equations of motion can be 

( 2 . 1 ) 


1 

It should be noted that in Eq . (2.1), M is constant in the case of 
laminar flow. For the sake of generality* Eqs. (2.1) and (2.2) can be 
nond.iinensionalized by introducing the following variables; 


written in tensor notation as 


Momentun ; 


Dt 


IP . 

3ic. 

1 


3X. 

J 


3u^ a u . 

. 3X. 3X. 

L J 1. 


Continuity; 


3U. 

1 


= 0 
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(2.3) 


A substitution of these variables in Eqs. (2.1) and (2.2) results in 


DU. 

1 


Dt 

3U. 

X 



3P . _3__ 

3X. 3X. 

J- J 


a 


3U. 3U. 

‘‘3X, 3X.^ 

J 1 


(2.4) 


(2.5) 


where a » p /Re (2.6) 

Equations (2.4) and (2.5) can be written in two-dimensional cartesian 
or cylindrical coordinates as 


3U 


3 U 


3U 


3t'^^3r‘^^3Z 


IZ L. 

3Z ^ 3Z 


r 


3Z’' 


h 


a r 


fUl_ 

'“3Z 3r'' 


(2.7) 


3 V 


3 V 


3 V 


3 t 3 r ^ 3 Z 


.3V 3U^ 2x6 /3V V> 


3P ^ 3 / 3V>> 

3 r ^ 3r 3r^* 

2x6 /3V 


3Z 


( 2 . 8 ) 




3U 

3Z 


= 0 


(2.9) 


Equations (2.7) to (2.9) are applicable to plane two-dimensional flows 
if 6 * 0 and to axisymmetry flows if 6 =1. 

2.2 Stream Function and Vorticity Formulations 
There are now three equations, Eqs. (2.7) - (2.9), and three 
unknowns U, V, and P. Introducing the definitions of the stream 
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funcfion and vorticity it is possible to reduce the number of equa- 
tions to tvo and this eliminates the pressure. The definitions of the 
stream function and vorticity are: 


3 V 3U 



where (u is the vorticity and i}/ is the stream function, 
A combination of Eqs . (2.7) through (2.12) results in 


( 2 . 10 ) 

( 2 . 11 ) 


( 2 . 12 ) 



32i|> 6 dijj 6 


= 0 


(2.13) 

(2.14) 


Equation (2.13) is the non-conservative form of the vorticity equation. 
The conservative fom can be obtained by multiplying the continuity 
equation by w and adding it to Eq. (2.13) as 
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The preceding equations are nonlinear parabolic-elliptic partial dif- 
ferential equations which can be solved numerically along with appro- 
priate boundary and initial conditions. 

2.3 Coordinate Transformations 
The governing Eqs. (2.13) through (2.15) are expressed in 
physical coordinates. For many problems, the boundaries may be quite 
irregular (Fig. 2.1a). This requires special consideration for the 
application of boundary conditions such as interpolation or some kind 
of approximations. Also a local grid resolution is required in most 
practical problems, which makes it extremely difficult to solve the 
governing equations in the physical coordinates. Therefore, it becomes 
advantageous to transform Eqs. (2.13) through (2.15) into new 
computational coordinates. The computational domain is an idealized 
rectangular coordinate system where a uniform grid is speci- 
fied (Fig. 2.1b). In other words, this transformation maps the z, r 
domain of the physical coordinates (Fig. 2.1a) into the ?, n domain of 
the new computational coordinates (Fig. 2.1b). This transformation, 
however, adds considerable complexity to the equations of motion. The 
following chain rules are used in the transformation process: 


f 

z 


z 5 


+ n 

z 


f 

n 


f 

r 



+ h 

r 


f 

n 


fp + 2 ? n fp_ + n f + f „ 
r rr 5 r r Cn rr n r nn 


rr 


f = ^ * 2 K n f,- f *h^f 

zz z zz K z z ' zz Ti * z nn 


(2.16) 

(2.17) 

(2.18) 
(2.19) 
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f 

zr 


-5 




* ^ zr h * r ^ ^ J ^Pr> * f + 

^ * r r z 5n zr n 


n n f 
z r nn 


( 2 . 20 ) 


Equations (2.13) through (2.15) now can be written in terms of the 
transformed variables as 


“ t ■ 4 “55 * “i “5 ♦ =1 “5n ♦ “i + % . 0 ^^ F; 

0 - 4 ♦55 . . P2 » 14 ^ ^ 


( 2 . 21 ) 

( 2 . 22 ) 


For the non-conservative form of the vorticity equation, Bq. (2.21), 
the coefficients are defined as follows; 

= a + ^2^) 

^ *“^^zz ^^rr * ^ J “ [u-2 ( 5 ^ 

-{ V-2 € a + q a )]^ 
r 5 r n ■' r 

q = 2X S ^ T1 ^ ^ n p 

^ zz rr /r) ^ J - [ U-2 (? , a + n a ) ]n 

^ 2 ^ z n'^-' z 

T. V-2 (S a + n a )]n 
r ? r n ■' r 

B a a (n^ + p2 ) 
z r 
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Fi - (fiw/r) [v + +’^r“n^ “ ~ 


55 


+(5 "-5 ) a. ■♦■ 2 (5 »i - 5 n 

2 z rr 5 z z r r 5n zz rr ti 


+(n2 ~ ) a^J [c V. + n + 5 u.+ n u ] 

z r nn z 5 z T) r 5 r n 

+2[? ? oc__ + C a + (c n +5 n ) a 

zrC^ zr5 zr rz 5tt 

+ n a +n n a J [c v + n v -C u - n u ] 

zr ri z r rin '■ r ^ rn z5 zn-' 


For laminar flow these coefficients will be reduced to the following: 

h = « (52^ * i\) 

Q. ^ 2^(5 *^r%> 

El +’^rr " ’^r " 

B = a (n^ + ) 

^ ' z r' 

FI = («<^/r)(V - ct/r) 

For the conservative form of the vorticity equation, Eq. (2.15), the 
corresponding coefficients are 


^ , a (^2^ + C2^) 

91 = “ 2Z + ^ rr + ^ ^ 

0 

G*2a^.n +5 n) 

^ 2 z r r 


r*': 
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^. -« +n2p 


„ _ rawfi . - 3(Vo) . _ 3(Vo) . , 9 (Ui») . _ 9(Uoj)i 

^1 " 1 — ■* 

t2 *^9? 9n 9? 9n 


The coefficients for the stream function eK^uatlon, Eq. (2.22), are 
defined as follows: 

% -5^, («/r)5r 

Co » 2(? n + ? n ) 

^ z z r r 

E2 - 

< z r 

6 

F 2 =» r w 

Equations (2.11) and (2.12) are written in the transfomed coordi- 
nates as 

U » (1/r)^ ^ r r ’^n ^ 

V = -(1/r/ „) 

The next step is to establish relations between the physical and compu- 
tational coordinates such that 
z = z(5 , T1 ) 

r = r(5 , h ) 


(2.23) 

(2.24) 


(2.25) 
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The above relations should be unique « single- valued and have continuous 
derivatives. This will be true if the determinant of the Jacobian 
matrix of Eq. (2.25) exists and is nonzero. The Jacobian matrix of the 
transformation is expressed as 


J " 





(2.26) 


The inverse transformation of Eq. (2.25) and its Jacobian matrix can be 
written as 

n n 

[ z ■ r J 

Equations (2:, 26) and (2.28) are related by 

j»[j*]-l (2.29) 


(2.27) 

(2.28) 


The relations between derivatives of the physical coordinates can be 
deduced from Eq. (2.29) as 


«r - -\T\ |J] 

where jj | is the Jacobian determinant. 


(2.30) 
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The relations between the second derivat,ives can be obtained by using 
Eqs. (2.16), (2.17) and (2.30) as 


zz 


rr 


rz 



(2.31a) 

"Cn 

(2.31b) 


(2.31c) 

*"z'5n 

(2.31d) 

+ C Zp- -% n Jp - j ) / | j j 
r r 5 r n ' * 

(2.31e) 

+ ? Tpp + n n J + ? n Jp)/I J 1 
r z r n r z 5 ' ‘ 

(2.31f) 


Equations (2.30) and (2.31) relate derivatives of physical coordinates 
to computational coordinates; these derivatives should exist. 
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Chapter 3 

COMPUTATIONAL PROCEDURES 

The governing equations, Eqs. (2. 13)”(2. 15) , are derived in 
Chap. 2 for incompressible flow in physical coordinates. These 
equations are transformed to new computational coordinates, as Eqs. 
(2.21) and (2.22). Aa mentioned before, these are full parabolic- 
elliptic partial differential equations, which are controlled by the 
boundary conditions for all variables along a surface which encloses 
the domain of interest. In case of turbulent flow, the equations are 
completed by supplying some auxiliary transport property relations. 
Because of the complexity of the equations there are only few 
analytical solutions. However, the equations can be solved using 
numerical techniques such as finite difference method or finite element 
method. There is also the method of lines (>£}L). In the present 
study, the stream-function equation is solved by the successive 
overrelaxation (SOR) method and the vorticity equation is solved by the 
method of lines, 

3.1 Stream Function Equatiif'^ 

The stream- function equation, Eq, (2.22), is a two-dimensional 
elliptic partial differential equations. Discretization of Eq . (2.22) 
yields a system of linear equations which should be solved at each 
iteration. There are several methods available to solve this equation 
some of which are listed below: 

1. Direct methods 
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Richardsoa's method (point iteration) 

3. Liebman's method (Gauss -Seidel) 

4. Successive overrelaxation (extrapolated Liebman method) 

5. Alternating directions implicit methods 

6. Hopscotch methods 

Since only the right hand side of this linear system, Eq. 

(2.22) , chaos<i8 at each iteration, it can be solved by the inversion 
method. This method sometimes is called the direct method. This is an 
efficient method for small systems, but the round off error destroys 
the accuracy of solution for large systems. 

The Richardson's method is also known as the Gauss iteration or 
point iteration method. As mentioned before, discretization of Eq. 

(2.22) yields a system of linear equations which can also be solved 
iteratively. The computation is based on values computed in the 
previous iteration, i.e., 

- f ['l»(N,k)j (3.1) 

where k is the iteration number. This iterative computation is carried 
out until some convergence criteria is satisfied. 

Liebman's method is like Richardson's method except it uses all 
updated -values at each iteration step, i.e.. 


= f ['P (N,k,k+1)] 


(3.2) 


I 




For better accuracy and efficiency^ the results from Liebman's method 
may be extrapolated as 


- ij;*' + u) f[ij/ (N,k,k+1)] (3.3) 

N 

where cog > 0 

In the preceding equation, is a relaxation parameter. If the value 

of is between one and two (l«ii^<2), it represents overrelaxation, 

and ifu> is less than one (o<u) <1), it represents und^irrelaxation. 

0 0 

Use of the over-relaxation tnethod usually is not recommended for those 
equations which contain strong source terms. It has been found 
experimentally [12] that is inversely proportional to Reynolds 

number. Therefore, underrelaxation is required as the Reynolds nimber 
is increased. 

Extensive discussions on alternating directions implicit (ADI) 
and hopscotch methods are available in [11-13] and are not discussed 
any further here. 


3.2 Successive Overrelaxation Formulation 
In this study the successive overrelaxation (SOR) technique is 
formulated in such a way that a one-dimensional array is used to 
compute and store each variable. A typical node arrangement for SOR is 
shown in Fig. 3.1. By using second order differencing, Eq . (2.22) is 
expressed as 


. r , k , k+1 , k+l-\ Bn r , k . k+li 

+ ['ll ^ 1 

7^‘-^N+IN+1 ^N+IN-1 ^N"IN-1 ^N-IN+1^ 


k+l T „ « 

* 't’K-INl* ^2“° 


^ _ r . h k+1 1 _ r k - k+l 

^ ^ L'J'N+IN "’^'n-InJ * ^ [’<'N+IN ■ 2 'f'N 


(3.4) 





ORIGINAL PAGE B 
OF POOR QUALfT^ 


Note that the coefficients are evaluated at the Nth node point. 
Rearrangement of Eq. (3.4) results in 

^ Residual (k, k+1, N) 

’^N “’^'n 2 (Ag + Ez) ^ (3.5) 

where the Residual (k, k-t-l, N) is defined as 


(2Ag + %) 2^2-82 

Residual (k, k+1, N) ^ ^ 

^‘•^N+IN+1 ^N+IN-1 ^N-IN-1 ^N-IN+1-' 

^ If ^ If If 

* r-*LiN " — <-m - ♦ Ea) 


(3.6) 


Equation (3.5) is Liebman's method, which can be improved by appropri-* 
ate extrapolation 


0) 

( k+1 _ ( k ° Residual (k,k+l,N) (3.7) 

^N " ^N 2(i^ + EgJ 


There is no analytical way to find the optimal relaxation factor C^^) 
for this case. However, one may use the optimal relaxation factor for 
the Possion's equation with Dirichlet boundary conditions [12]. This 
is given by 



(3.8a) 
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where , 




+ COS 


(J ! — ) 


(3.8b) 


In Eq. (3.8fj), IN and JN are the number of grid points in the axial and 
the vertical (radial) directions respectively. 

3.3 Application of the Method of Lines to the Vorticity Equation 
Ifi the vorticity equation, the spatial derivatives are replaced 
by a corresponding set of second order difference equations. Also, the 
vorticity is considered to be continuous in time. Ihis gives 

d u) 

"dF” “ ^ ('^N+l ' ^N ’*'n- 1^ ®l('<^N+l ■■ '^N-l)'^^ 

^’*'n+IN+1" ’*'n+IN-1 ■^'^N-IN-1 " ’^M-IN+1^/'^ 

('j' jj+IN" 'I' n-In) (^’^N+IN ” ’I^N-In) * (3.9) 


Equation (3.9) is a set of coupled ordinary differential equations, 
which should be integrated simultaneously. 

The method of lines is used for solving the system of partial 
differential equations. As mentioned in the introduction, the method 
was developed and used originally in the Soviet Union. Liskovet's 
review article [1] is a review of 154 papers which date up to the mid- 
sixties. This review has considered linear partial differential 
equations' of elliptic, parabolic, and hyperbolic types. 

The method of lines basically changes a partial differential e- 
quation to a set of ordinary differential equation which can be 
integrated numerically. In the case of time dependent partial 
differential equations, the procedure is to discretize the spatial 
variable and treat the temporal variable as continuous in time. This 
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semi-discretization results in a set of ordinary differential equa- 
tions, which can be integrated along the lines of time. 

3.3.1 Advantages and Disadvantages of the Method of Lines 

Finite difference and finite element techniques are well-develop- 
ed and in general they are more attractive in terms of efficiency. 
Solutions of partial differential equations using these two methods 
lead to a system of linear algebraic equations which can be solved 
by direct or iterative methods. However, the nonlinear teras in the 
equations must be linearized. In the method of lines, there is no need 
for linearization. Other advantages are as follows: (1) possible 

establishment of convergence and stability criteria for a wide class of 
problems [3, 14] , (2) accurate solutions by using higher order 
approximations for the derivatives and with less computational costs in 
comparison to the finite difference method, and (3) more efficient due 
to a better time-step control and easier implementation even for non- 
linear systems. 

There are certain limitations in the method of lines such as the 
number of lines. Jones and et al. [5] estimated that error is equal to 
exp [ 4l^^ A b] , where N is the number of lines and b is the characteris- 
tic length. This means that using a large number of lines (for better 
resolution) may bring significant instability. The former conclusion 
is valid for elliptic systems. 

3.3.2 Convergence, Accuracy, and Stability 

Convergence exists when the solution approaches the solution of 
original continuum differential equations as step size or grid size 
approaches zero. But, instability occurs when round off error or any 
other computational errors become unbounded. There are extensive 


discussions about convergence and stability by Jones et al. [5], and 
Loeb and Schiesser [10]. Their results show that using higher order 
finite difference approximations for derivatives improves the accuracy 
and the stability, and reduces computer cost. They have also shown 
that using more grid points improves the accuracy, but makes the 
solution become less stable. Stability may be related to the stiffness 
of the equations. This effect is investigated in the next chapter. 


Chapter 4 

STIFFNESS ANALYSIS 


In engineering problems , stiffness may arise due to the physics 
of the problem or be introduced to the problem because of the type of 
approach applied for obtaining solutions. An example is a problem in~ 
volving chemical reactions where time scales span from 10"® to 10® 
seconds, simultaneously. This system is referred to as a stiff system 
when the processes are coupled and when all time scales must be resolv 
ed. In mathematical terms, when the eigenvalues of a system of 
differential equations have a large variation, the system is referred 
to as a stiff system. For example, consider the following equation: 


y *Ay, where A “ 


-1 

0 

0 


49 

-50 

1150 


0 

0 

-1200 


(4.1) 


This equation has a solution of the form 


Yi (x) 
^2 (x) 
Yg (x) 


-X “50x 

- e + e 

-50x 

= e 

“50x -1200x 

= e + e 


The eigenvalues of this system vary from -1200 to -1. The degree of 
stiffness is measured by the stiffness ratio which is the ratio of the 
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largest eigenvalue to the smallest eigenvalue of the system. For 
example, the stiffness ratio for Eq. (4.1) is 1200. Stiffness may be a 
local problem in which case the equations are stiff in some regions and 
nonstiff in other regions. Numerical integration of Eq. (4.1) requires 
a specified integration step size h which is determined only by the 
components of the eigenvalues, and the stability region of the inte- 
gration scheme. Ihe stability of numerical integration of Eq . (4.1) is 
governed by the maximum absolute value of X^. For example, using 
Euler's method for the integration, it is necessary that {Max (X^)}*h<2 
[15]. This implies that the maximum stable step size is 1/600, meaning 
that 600 integration steps are required to reach x = 1. 

There are only a limited nauber of numerical methods available to 
solve ordinary differential equations which utilize the stiffness char- 
acteristic [9, 14-17]. A popular method is the Gear's method [9, 14, 
17]. Subroutine VOADAM [9] has been used to solve Eq . (3.9), which 
has the options for stiff and nonstiff systems. Stiff and nonstiff 
options differ in storage and computer time. The computer memory 
required for a stiff solver is of order (N nunber of equations), 
whereas for a nonstiff solver it is of order N. However, the stiff 
option generally requires much less computer time. As an example, 
Madsen and Sincovec [6] observed a 400% computer time saving by using 
a stiff integrator. 

4.1 Effect of Grid Concentration on Stiffness 

In many physical problems, one desires to have an arbitrary grid 
distribution with concentration at any location in the physical domain. 
In particular in flow field computation, concentrations are required 
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to capture rapid changes in flow conditions such aa boundary layer, 
shock, and separation. Concentrations affect the stiffness 
characteristics of a system, and it is essential to investigate the 
effects of grid concentration or stiffness. For this purpose, a simple 
one-dimensional heat equation is considered. Ihe equation is given 
by: 


3 T 3^T 


0 < X < 1 


(4.2) 


Equation (4.2) is transformed from the physical coordinate (x) into a 
computational coordinate (n ) 


3 T 
3 t 



+ n 


XX 




(4.3) 


where n“h(x), o<ri<l 


(4.4) 


In Eq. (4.3), the spatial derivatives are discretized using a second- 
order finite-difference approximation and T is assimed to be continuous 
in time. This results in 



a. T. , 
1 1+1 


+ b. T. + c. T. , 
11 1 1-1 


where 


(4.5) 


(4 . 6a) 


b. = 2(n /An)^ (4.6 b) 

c . “ (n /2&n) (4.6c) 

XX 
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Equation (4.5) represents a set of ordinary differential equations 
which can be written in matrix form as 


dt 


[A T 


(4.7) 


where [ a] “ 


bi Cl 

h < ^ 


I 


a b c 

N-1 N-1 ^N-1 


Eigenvalues of matrix [ a] are determined in order to analyze the stiff- 
ness characteristics. Eigenvalues of matrix [a] depend on boundary 
conditions and the relation between physical and computational coordi- 
nates. 


Consider the following relation between x and n, and the boundary 
conditions: 


X =» b( e^ -l) 

(4.8a) 

K 


b = ( e -l) 

(4.8b) 

T(0) - — ^ » 0 

(4.8c) 


In the above equation, K is a stretching factor which determines the 
degree and location of the grid concentration. A large positive value 
of K results in a high grid concentration near x = 0, whereas a large 
negative value of k results in a high grid concentration near the x = 1 
boundary. A substitution of Eq. (4.8) into Eq. (4.6) results in 
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a. 

1 


-2»i ^ -;>wi ^ 

e __ e 

(&n Kb)^ 2Kh?An 


-2»1 . 

j,. .Z2e 

^ (An Kb)^ 


c. 

1 


-2»1 , -2W1 . 

s— t 

(AriKb)2 2Kl^An 


(4.9a) 


(4.9b) 


(4.9c) 


Eigenvalues of matrix [a] are found numerically using subroutine REQR 
[18], and the stiffness ratios are computed for different values of 
stretching factors. Figure 4.1 shows a plot of the stretching factor 
versus the stiffness ratio. The plot indicates that stiffness in- 
creases with the magnitude of the stretching factor. Also, the slope 
of the stiffness characteristic curve is higher near the fixed boundary 


(x"0) than near the derivative boundary (x=l). 

4.2 Effects of Differencing Scheme on Stiffness 
In the finite difference approach, convergence and stability of a 
solution depend on the differencing scheme used. Similarly, the 
differencing scheme affects the stiffness of the resulting ordinary 
differential equation system in the method of lines. To investigate 
the effects of differencing schemes on stiffness, consider the 
linearized one-dimensional advection-dif fusion fluid flow equation: 


do) 

a t 


3bl 



1 3^0) 

Re al?" 


(4.10) 
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(i}(0,t) “ 0 

w(l,t) - I 


ORSGK^At - - _ 

OF POOR QUAUn 


( 441 ) 


where 0 < x < 1 


In Eq. (440), the spatial derivatives can be discretized by using the 
backward,, central, or forward differencing scheme. Thus, Eq . (4.11) 
can be expressed as 


di) i 
dt 


a w i-1 + b t») i + c OJ i-fl 


(4.12a) 


This equation is a set of ordinary differential equations which can be 
expressed also in matrix form as 


do 

dt 




where 


[a] 


b c 

a b 


c 

b 

a 



(4.12b) 


(4.13) 


The values of elements a, b, and c of matrix [a] are functions of the 
nunber of grid points (N) and the Reynolds number (Re) .For forward 
differencing, they are 

a = I?/Re , b » N-2I^/Re, c = N^^Re-N 

For backward differencing, they are 

a « /Re + N , b » -N-2I^/Re , c » N^/Re 
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For central differencing, they are 


a » /Re + N/2, b « -21^ /Re , c “ N'VRe-N/2 

The eigenvalues of matrix [a] can be computed analytically, and they 
are 

An - b + 2 VaT COS (4,14) 

A systan of ordinary differential equations (e.g., Eq. 4.12) is 
unstable if the real parts of its eigenvalues are positive. Iherefore, 
backward and central differencing are inherently stable, whereas 
forward differencing is conditionally stable, provided Re < 2 N. It is 
shown later that forward differencing is not useful for solution of Eq. 
(4.10). 

Stiffness ratios for these differencing schemes are computed and 
are plotted on Fig. 4.2. The plots show that the equations become 
stiff as the number of grid points is increased (curves 1-3), and the 
stiffness ratios do not depend on the differencing scheme at low 
Reynolds number (curves 4 and 5) 

A specified integration step size is required for the integration 

of ordinary differential equations. Maximum allowable step size 

depends on the eigenvalues and the applied integration techniques. The 

step size should be selected such that Ah is located in the stability 

region (Fig. 4.3). For instance, the stability region is shown in Fig. 

4.3 for the Adams -Bash forth method (first order) where h is the step 

size and A is the eigenvalue. Generally, the integration step size 

is inversely proportional to the distance between the location of the 

eigenvalue A and the stability region. For example, if A , A , and 

A D 





Fig. 4.3 Stability region for Adams-Bashforth (first order). 



X are the largest eigenvalues of three different systems of ordinary 
differential equations (Fig. 4.3), system A requires the smallest step 
size, because X has the farthest distance from Che stability region. 

A 

The maximum allowable step size can be computed analytically for Che 
Adams-Bashforth technique (first order). Consider eigenvalue X^ * 
RE(^ ) + i IM(X)], where RE is the real part and IM is the imaginary 
part, the step size should be selected such that the point A in the 
complex plane can be moved to the stability region. This requires 

h^ [Modulus of X *] = [Modulus of X 1 (4.15) 

where from Fig. 4.3 Modulus of X^ = AO, Modulus of Xa = aO. 

After some algebraic manipulation, Che maximum step size is found to be 

h ■ — — — - (4.16) 

^ [Modulus of X 

Using this equation, the step size can be computed for the different 

differencing schemes. In case of forward differencing, the step size 

is positive if Re > 2N, and the system is unstable. Therefore, it is 

impossible to solve Eq . (4.11) using forward differencing along with 

the Adams-Bashforth (first order) integration technique. Similarly, 

step size for central differencing is 

h„ “ 1/Re. (4.17) 

Max 

The step size is inversely proportional to Reynolds number for Re > 2N. 

In case of backward differencing, the step size is 

h “ 2/ (Modulus of X ) (4.18) 

Max 

Figure 4.4 shows variations of step size with Reynolds number (Re) and 
the number of grid points (N) . The figure shows that the step size de- 
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creases as the number of grid points increases (curves 1 and 2). Ihis 
is compatible with the finite difference techniques [12]. On the same 
figure, step size is plotted versus Reynolds number. For backward 
differencing, the step size increases with increasing of the Reynolds 
mmbec. The converse is true for central differencing. 

For inviscid flow (Re-»- ®), forward differencing is inherently 
unstable. For central differencing, step size is zero which makes it 
impossible to integrate the equations. The step size for backward dif- 
ferencing is 

h = 2/N (4.19) 

The step size is inversely proportional to the number of grid points; 
this result is similar to the result of viscous flow. 


Chapter 5 


BOUNDARY AND INITIAL CONDITIONS 

For a particular flow problem, selection of proper boundary 
conditions depends upon the nature of the flow and the computational 
procedures employed. The application of certain conditions may cause 
nimerical instability in the solution even though the flow is physic- 
ally stable. 

Equations (2.21) and (2.22) are parabolic-elliptic partial dif- 
ferential equations. The dependent variables in these equations should 
be defined by some relations along the boundaries. There are three 
general types of boundary conditions for a dependent variable and they 
can be stated as follows: 

1. Specifying values of the dependent variables at the boundar- 
ies . 

2. Specifying first or higher derivatives of the dependent vari- 
ables at the boundaries . 

3. Specifying algebraic relations which relate dependent vari- 
ables . to independent variables or to their first or higher 
order derivatives!. 

Two points should be noted in choosing the above conditions. 
First, the second condition cannot be applied at all boundaries because 
uniqueness considerations. Second, a combination of the above 
conditions can be applied to various parts of one boundary. 
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In the selection of boundary conditions, four important factors 
that should be considered are: convergence, stability, computer time, 
and above all, the physical justification. Extensive discussions of 


these factor are given by Roache [12]. 

For the confined flow geometry, Fig. 5.1a, there are four types 
of boundary conditions inlet: outlet, wall, and synnetry. Whereas, in 

the driven cavity. Fig. 5.1b, there are two types of boundary condi- 
tions: wall and moving wall. Also, in each case, it is necessary to 

specify some initial conditions. 

5.1 Inlet Condition 

For problems involving duct flow, the inlet conditions are 
usually specified. Throughout this study, the inlet conditions are 
fixed. In general inlet velocity profile for duct flow is given by 

U(r) « 1 - AI*r^^^, (5.1) 

V(r) = 0.0 (5.2) 

The stream function is computed from Eqs. (2.11) .and (2.12) by 


6+1 NI+6+1 

^ “ 5TI TTi-hS +1 


(5.3) 


The vorticity is computed similarly from Eq. (2.10) as 

NI-’l 

0) = AI • NI • r 


(5.4) 


In the preceding equations, NI and AI are the inlet distribution para- 
meters and they can be selected to produce desired inlet conditions. 
For example, AI = 0 corresponds to a uniform distribution of velocity 
across the inlet; AI = NI “ 1 corresponds to a linear distribution of 
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inlet profile; AI * 1, NI “ 2 corresponds to a parabolic inlet profile. 
In this study, a uniform inlet profile is selected for flows in pipe 
and parallel ducts, and a parabolic profile for flows in a curved~wall 
diffuser. 

5.2 Outlet Condition 

For flow between parallel plates and in pipes, the most realistic 
outlet condition would be no profile changes in the flow direction far 
from the entrance point; provided there is no change in the wall 
conditions. Ihis makes it impossible to use realistic conditions for 
such cases. Roache [12] points out that numerical experiences show 
that catastrophic instability may be propagated upstream from the 
application of improper outlet conditions and this may destroy the 
solution completely. 

The Reynolds number is an important criteria for selection of an 
outlet condition. For high Reynolds number flows, the derivatives with 
respect to flow direction are generally small compared to the deriva- 
tives with respect to normal direction. Therefore, the governing equa- 
tions tend to be parabolic in nature (bound ary -layer equations) . In 
this case, the outlet condition has little effect on the solution. 

But, low Reynolds number flows require a physically well justifiable 
outlet condition. Generally, the outlet condition can be described 
to be fixed or to be time dependent. A fixed outlet condition is the 
easiest to apply from a computational viewpoint. For a fixed outlet 
condition, the velocity profile is given by 

U(r) = U (1-AIr**^) . 

max 


(5.5) 

(5.6) 
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The streata function and vorticity can be computed in the same way aa 
described in Sec, 5.1. They are 


6 -fl 


0) (r) 


U I i-z AI » , 

max '■ 6 +1 NI-*^ +1 


NI+6 +1 
' ■] 


U AI 
max 


NI 


NI-1 


(5.7) 

(5.8) 


In the absence of suction and/or blowing, the stream function should be 
constant along the wall. This is given by 

(5.3) ' (5.7) <5.9) 

A combination of equations (5.7) and (5.9) results in 

U = 1/(S+1) - AI/(NI + S+L) (5.10) 

max ^+i//r N0-hS+l'( 

/(6+l)-A0(r^ J/(N0+6+1) 

where r^ is the outlet radius. 

A fixed outlet condition is not suitable for separated flows or 
any flows with a viscous wake [12]. Paris and VRiitaker [19] solved a 
confined flow using zero gradient conditions. This showed improvement 
in solution over .specified conditions. Later, Thoman and Szewczyk [20] 
used physically less restricted conditions for the stream function and 
zero gradient on the vorticity as 

31U „ (5.11) 

® 

Briley (21) took a different approach by considering the following re- 
lation 

(5.12) 
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A subacitution of Eq. (5.12) into Eq. (2.14) results in 

iM < 5 . 13 ) 

“ - a/r) [JT- -j^J 

Equations (3.12) and (5.13) represent two ways of conq>uting vorticity 
at an outlet. This technique is referred to as maltreatment of outlet 
condition by Roache [22]. Throughout this study, Eq. (5.11) is 
considered for the outlet condition. 

5.3 Wall Condition 

An impermeable wall assunption allows the stream function to have 
a fixed value along a wall regardless of its geometry. The stream 
function may be computed from the inlet condition, Eq. (5.3), 

AI (5.14) 

^ 6 +1 " III-hS +1 


The vorticity can be computed by applying Eq. (5.14) to Eq . 
(2.22) , i. e. , 


^ + Ej + F2 =0 


(5.15) 


.Also, zero velocity at the wall allows 


» 0 
n 


(5.16) 


Upon combining Eqs. (5.14) and (5.15), one obtains 

+ r\2 


Hi 


n*- n- 

z r 




(5.17) 


nn 
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Equation (5.17) describes variations of vorticity along a horizontal 
wall. K similar expression can be derived for a vertical wall in the 
case of a driven cavity, i.e., 

r2 + r2 

<5.18) 

o 

r 

Equation (3.17) is applied to the boundary A1 shown in Fig. 5.1a and to 
boundary B3 shown in Fig. 5.1b. Equation (5. .18) is applied to bound- 
aries B2 and B4 in Fig, 5.1b. 

5.4 Symmetry Condition 

Selection of a symmetry boundary condition depends on the nature 
of the flow. In this study, since there is no center body, the 
syntnetry condition along the centerline requires that 



(5.19) 


V » 0 (5.20) 

By comparing Eqs. (5.18) and (2.12), it can be concluded that the 
stream function has a fixed value along the line of symmetry which can 
be chosen to be zero. Similarly, it can be concluded that the 
vorticity is zero along the line of symmetry. 

5.5 Moving Wall Condition 

In Fig. 5.1b, the Bl boundary moves with a uniform velocity U=-l. 
This means that the boundary can be chosen as a streamline. Ihe value 
of the stream function can be selected to be zero to match with the 
value of the stream function at boundaries B2 and B4. Using uniform 
velocity at the boundary, Eq. (2.23) can be written as 

? 'I'j, + T) “ -1 (5.21) 

r 5 r n 
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By using the value of the stream function along the boundary (<p * 0), 
this equation reduces to 

tn ■ -1/n^ (5.22) 

Vorticity can be obtained from Eq, (2.22) by using the value of the 
strean function along the boundary and Eq. (5.22) | i. e., 

w » Dlz/n^ - ’1'^^ (5.23) 


The second derivative of the stream function can be approximated 
(second order) by 


'I' 


nn 


N+IN 


2 


N-IN 


(5.24) 


Also, the first derivative of the stream function can be approximated 
(third order) by 


(2|/ 


N+IN 


^N- 


th + 

^N“IN 


'^N-2IW^ 


-1/n. 


(5,25) 


The combination of Eqs. 


(5.24) 
N-IN ” 


and (5.25) results in 
N-2IN - ^ 


(5.26) 


Equation (5.26) is a second order accurate equation [23]. 

A combination of Eqs. (5.23) and (5.26) results in 

W - + Ez (6/n^ - 8 + 4 »j^_2IN ^ (5.27) 

A similar approach has been employed for evaluating the vorticity at 
the stationary wall boundaries. 
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5.6 Inilial Conditions 


In computational fluid dynamics, the initial conditions usually 
correspond to a real initial situation for a transient problem or, as a 
rough guess, for a steady state problem. In practice, initial condi- 
tions are obtained from experiments, empirical relations, or approxi- 
mate theories. The initial conditions used in the determination of the 
steady state solution should have no significance in the steady-state 
solution of incompressible flows [12], otherwise the solution is not 
unique. 

Generally, there are two kinds of initial conditions. In the 
first one, an impulse motion starts from the rest, and in the second 
kind, the flow has the same initial motion everywhere except the 
boundaries. In the present study, the entire flow-field is initially 
set equal to the inlet condition for the internal flows. For the case 
of driven cavity, the stream function and vorticity ar^ initially set 
equal to zero in the entire flow-field. 


Chapter 6 


GRID GENERATION TECHNIQUES 

It is highly desirable to solve partial differential equations in 
a rectangular box with uniform grid spacings [ 15 - 21 ]. Ihis is espe- 
cially true in fluid dynamics where the governing equations are com- 
plex. Ideally, a physical domain should be transformed to a computa- 
tional domain, where the physical boundaries map into the boundaries of 
a rectangle. This transformation has certain advantages. In the first 
place, the boundaries can be represented more accurately. Secondly, it 
allows better resolution in regions where rapid changes occur, such as 
boundary layers, shocks, and separated flows. Above all, computer 
codes which are applicable to any geometry can be written without the 
need of special procedures for the boundaries. The disadvantages are 
that the convergence, stability, and stiffness characteristics of the 
equations are affected. Also, the transformed equations are more 
complex then the original equations. The transfoirmed governing 

equations contain the rate of change of the computational coordi- 
» 

nates with respect to the physical coordinates. These derivatives are 
computed from the relations of the physical grid to the computational 
grid. There are, generally, three approaches for grid generations: 
classical technique (conformal mapping) f24J j differential methods 
[25], and algebraic methods [26]. 
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The classical technique requires the use o£ conq>lex functions to 
define the mapping which is extremely difficult for arbitrary geom-' 
etries. But, conformal mapping has the advantage of minimizing the 
number of terms in the transformed equations. The algebraic and 
differential methods for grid f^eneration are outlined in the following 
sections. 

6.1 Algebraic Method 

In the algebraic method an explicit functional relationship 
between the computational domain and the physical domain is determin- 
ed. The computational domain is rectangular and has a uniform grid 
distribution. The physical domain is defined by 

z ■ z (C , n ) (6.1) 

r = r (5 , n ) (6.2) 

A requirement of boundary-fitted cooidinate systems is that the bound- 
aries of physical domain map to the boundaries of the comp'utational 
domain. That is , 


Zl = z^ (5 , 0) = z^ 

(?) 

(6.3) 

11 = ri (C , 0) = ri 

(?) 

(6.4) 

li 

II 

(?) 

(6.5) 

i"2 “ T2 (? , 1) * 12 

(?) 

(6.6) 


Equations (6.1) and (6.2) can be rewritten using Eqs. (6.3) 
through (6.6) as 

z 3 z (5 , Ti ) » z (zj ,Z 2 ,n ) (6.7) 

r = r (5 , n ) = r (ri ,r 2 ) (6.8) 

The explicit forms of Eqs. (6.7) and (6.8) can be written as a simple 
parametric linear or cubic polynomials as [26] 
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Linear: 


z-a^ (?)n + Zi (5)(i-n) (6.9) 

r*i2 (C )n + ri(?)(l-J1) (6.10) 

Cubic: 

z ■ ^ (C ) fi (n) +s^ (5) ^ (n) + ^ f3<n) + 

dag (g ) fit (T1) 
on 


r-ri (?) fi (n) + 12 (?) ^ (n) + f3 (n) + 


dg> (?) % (n) (6.12) 

dn 


In Eqs. (6.11) and (6.12), the f's are the blending functions and are 
defined as 



- ^3 - 3h+l 

(6.13) 

h 

= + 3H 

(6.14) 

% 

= " 2t]^ -»T1 

(6.15) 


= 

(6.16) 


Equations (6.9) through (6,12) are called connecting functions. If f 
and n are the computation coordinate functions, then ? (?) and n(n) can 
be defined to produce a desired grid spacing distribution. Ihat is 

r-r(C) , 0<?<1 , Ocfci (6.17) 

tT =tT(ti) , CKri< 1 , O^n*; 1 (6.18) 


For example, contracting the physical grids towards one boundary can be 
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accomplished by using the following function 


Kn 



(6.19) 


e 

where K is a stretching parameter which determines the location and 
degree of concentration. 

This method (i.e., the algebraic method) requires relatively few 
computations. 


6.2 Differential method 

If the computational coordinates 5(z,r) and r|(z,r) are harmonic 
then the Jacobian does not vanish [27]. This allows a computation of ^ 
and n from an elliptic system (Laplace's equation), i.e., 

V2? - 0 (6.20) 

72n = 0 (6.21) 


The spacing of the coordinate lines can be controlled by adding 
inhomogeneous term to the right sides of Eqs. (6.20) and (6.21) as 

= P (? ,n) (6.21a) 

= Q (5 , n) (6.21b) 


In Eqs. (6.21a) and (6.21b), D and ? are known, and the unknowns are z 
and r. In order to be able to solve for z and r, the dependent 
variables of Eqs. (6.20) and (6.21) should be interchanged. Thus, 




(6.23) 
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where 


a 

0 

Y 

J 




5 n 


+ r r_ 

n ? 


+ r2 


Z_ r - z r_ 

? n n 5 


(6 . 24a) 
(6.24b) 
(6.24c) 
(6.24d) 


There are four boundary conditions for Eqs . (6.22) and (6.23). Tha 
quantities P and Q are the forcing functions which are selected to have 
a desired grid distribution and orthoganality at the boundaries [27 j. 
The disadvantage of this method is that it is inefficient compared to 
the algebraic method. 
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Chapter 7 


PHYSICAL APPLICATIONS 


This study investigates the validity and viability of the method 
of lines for solutions of the Navier-Stokes equations with arbitrary 
boundary geometry and grid spacing. For incompressible viscous flows, 
the pertinent equations are derived in two-dimensional and axisynmetric 
coordinate systems, Eqs. (2.13) -(2.15). These equations are trans- 
formed from the physical domain with arbitrary boundaries and grid 
spacings to a rectangular computational domain with uniform grid 
spacings. The resulting equations, Eqs. (2. 21) and (2.22), are solved 
numerically. The method of lines is used to solve the vorticity equa- 
tion, Eq . (2.22), and the successive overrelaxation technique is used 
for solving the stream function equation, Eq. (2.21). Boundary 
conditions at a solid stationary wall are no-slip conditions which 
corresponds to a fixed stream function and requires that the vorticity 
be computed from the stream function by Eqs. (5.17) and (5.26). For 
the internal flow computation, the stream function and vorticity are 
equal to zero along the line of symmetry. The inflow boundary 
condition is kept fixed throughout this sjtudy.^ The outflow condition 
is used to enforce no-change in the stream function and vorticity with 
respect to flow direction, Eq. (5.11). In the case of the driven 
cavity, the vorticity at the moving wall is computed from Eq . (5.28). 
For internal flow, the initial conditions are the inflow conditions for 
the entire domain, and no flow in the case of the driven cavity. Grid 


50 


distributions are generated by the algebraic method which makes it 
possible to handle any boundary geometries [26]. Specific solutions 
are obtained for internal flows (pipe, parallel plate and curved wall 
diffusers) and a driven cavity. These problems are discussed briefly 
in this chapter. 

7.1 Internal Flows 

The entrance regions of ducts are used usually to compare numeri- 
cal procedures, because of the availability of analytical solutions. 
Most of these solutions are based on parabolic forms of Navier-Stokes 
equations which are called bound ary -layer equations. Assumption of the 
boundary layer (neglecting the axial diffusion) makes it possible to 
obtain analytical solutions [28-31]. 

Vrentas et al. [32] investigated the effects of axial diffusion 
of vorticity on flow development in circular ducts. It was concluded 
that at very low Reynolds numbers (Re<20) , the axial diffusion causes 
the velocity development to be spread out in the downstream region, and 
this results in a larger entrance length than that predicted by the 
bound ary- layer analysis. It was also concluded that the converse is 
true for Re>20. The flow in parallel ducts may, in general, be divided 
into three regions (Fig. 7.1). Region I, is a region of potential 
flow. This flow interacts with region II which is a region of 
bound ary- layer flow. The third region (region III) is the region of 
fully -developed flow. Using these concepts, Schlichting [28] obtained 
an approximate solution by using Blasius' solution for the inlet 
regions between parallel ducts and a perturbation of the f ully-develop- 
ed solution in the region far from the entrance. The solution for the 
intermediate regions was obtained by patching together the upstream and 
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downsCream solutions. A similar approach was implemented by Atkins and 
Goldstein [33] for a circular tube. Later, Van Dyke [31] improved 
Schlichting' s solution by using more terms in the perturbation expansion. 
Another common method is the linearization of the convective term in the 
boundary-layer equations; this allows the velocity to be continuous along 
the axial coordinate. Using this concept, Sparrow et al. [30] developed 
solutions for circular tubes and parallel plates using bound ary -layer 
equations. In the present study, results obtained from the method of lines 
are compared with results of Sparrow et al. [30]. 

In the case of curved-wall diffusers, the outer wall geometry is 
based on the following equation: 

r - { 1 + (r2 - 1) ^ [ 1 + a(l - ~)]}H(X-L) + X H(L-X) 


where (7.1) 

V 2 = L tan 0 +1 

Figure 7.2 shows a typical curved-wall diffuser geometry. It is a 
bell-type diffuser for a>0 and a trumpet type diffuser for a<0. The 
value of a has been set equal to unity for the present study. By using 
the method of lines, solutions are obtained for different diffuser 
angles (6 ) and for a Reynolds number of 200. 

7.2 Driven Cavity 

The flow in a driven cavity has been used as a test case for 
evaluating the numerical procedures for solving the Navier-Stokes equa- 
tions. It serves as an example of closed stream lines which represent 
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interesting characteristics such as vortex development in the core and 
boundary- layer development on the walls of the cavity. 

Kawaguti [34] performed one of the earliest calculations for flow 
in a driven cavity and satisfactory results were obtained with Reynolds 
number up to 64. Later, Burggraf [35] obtained solutions for Reynolds 
numbers up to 400. An underrelaxation method was used to solve the 
governing equations, and a relaxation factor of 0.4 was used for a 
Reynolds number of 400. Burggraf was not successful in obtaining 
solutions for higher Reynolds numbers. Ihis may be due partially to 
the form of the governing equations (non-conservative form) used. 
Bozeman [36] solved the same problem using the conservative form of the 
Navier-Stokes equatious. Solutions were obtained successfully up to 
Reynolds numbers ct tOOO. Bozeman used a strongly implicit procedure 
which is unconditionally stable. Smith and Kidd [11] obtained 
solutions up to a Reynolds number of 5000 by solving the vorticity 
equation using alternating-direction implicit (ADI) methods and the 
stream function equation by the Buneman direct method (BDM) . Using 
the method of lines, Kurtz et al. [8] were able to obtain solutions up 
to Reynolds numbers of 50,000. In this study, solutions are obtained 
for three Reynolds numbers (100, 1000, and 10000) with two different 
16x16 grid distributions. 
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Chapter 8 


ElESULTS AND DISCUSSION 


By using the formulation and numerical procedure discussed in the 
preceding chapters, solutions are obtained for flow' between parallel 
plates, circular pipes, curved-wall diffusers, and in a driven cavity. 
Non-conservative forms of the governing equations are used for the 
entire internal flow computations. 

Figure 8.1 shows the 20x20 grid distribution which is used for 
computing the flow field between parallel plates and in pipes. Ihe 
grids are concentrated near the walls and the entrance point where 
rapid changes occur. 

For a Heynolds mjonber of 200, Figs. 8.2 and 8.3 show velocity, 
stream function and vorticity distributions for flow between parallel 
plates and in a pipe, respectively. Velocity profiles are con^ared 
with those obtained by Sparrow et al. [30] , and are found to be in good 
agreement except in the vicinity of the entrance point. This discrep- 
ancy may be due to the fact that Sparrow's results ai:e based on the 
bound ary- layer assunptions. Ihere is a sharp drop in the stream 
function in the vicinity of the entrance point which indicates a rapid 
flow development in this region. Uiis sharp drop corresponds to a 
sharp change in vorticity which can be observed in Figs. 8.2c and 8.3c. 
Similar results are obtained for a Reynolds number of 1000 and these 
are illustrated in Figs. 8.4 and 8.5 for flow between parallel plates 
and in a pipe, respectively. 
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Fig. 8.3a Velocity distribution for pipe flow. Re = 200. 
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Fig, 8.4b Stream function contours for flow between parallel 
plates. Re = 1,000. 
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Fig. 8.4c Vorticity contours for flow between parallel plates. 
Re = 1,000. 
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Fig. 8.5a Velocity distribution for pipe flow, Re 









Figure 8.6 shows a 20x20 grid distribution for a five~degree 
curved-wall diffuser. Again, the grid is concentrated near the en- 
trance point and wall region for capturing the boundary layer develop- 
ment and separation. Figure 8.7 shows the velocity profiles, stream 
function and vorticity distributions for a Reynolds number of 200. The 
results show no indication of separation. The next case studied is a 
ten-degree curved-wall diffuser. The grid distribution for this case 
is shown in Fig. 8.8, and the results shown in Fig. 8.9 clearly 
indicate that a separation occurs in this case. 

Figure 8.10 shows a uniform 16x16 grid distribution for a driven 
cavity. Figure 8.11 shows the the stream function and vorticity 
distributions for Reynolds number of 100. These results are obtained 
using the non-conservative form of the Navier-Stokes equations with a 
uniform grid distribution. The results are in good agreejnent with the 
results obtained by other investigators [8, 11, 23, 35-37]. For a 
Reynolds number of 1000, the solution does not converge if the non- 
conservative fomn of the governing equation is used. For higher 
Reynolds numbers, therefore, the conservative form of the governing 
equations, Eqs. (2.15), are used to obtain further results. Figure 
8.12 shows the stream function and vorticity distributions for a 
Reynolds number of 1000. The rs'i^isults are in good agreement with Kurtz 
et al. [8], Smith et al. [11], Ghia et al. [23], Bozeman [36], and 
Schreiber and Keller [37]. Ghia and Schreiber have used over 10** grid 
points. Figure 8.13 shows the stream function and vorticity distribu- 
tions for Reynolds number of 10,000. The results are in good agreement 
with the results of Kurtz et al. [8], Ghia et al. [23] and Schreiber 
and Keller [37]. These results suggest that the high Reynolds number 
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Fig. 8.7b Vorticity contour for a bell-type diffuser 
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8.7c. Stream function contour for a bell-type diffuser 
a = 1, 0 = 5, L = 20, = 32, = 200. 
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Fig. 8.9c Vorticity contour for a bell-type diffuser 
0 = 5, L = 10, Lj. = 40, Re = 200. 
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Fig. 8.10 Uniform grid distribution for 


a driven cavity. 
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Fig. 8.13b Vorticity distribution for uniform grid 
distribution, Re = 10,000. 




flow conaists eaaentially of a single inviscid core of vorticity with 
viscous effects being confined to a thin shear layer near the wall 
boundaries. By using a nonuniform grid distribution, Fig. 8.14, 
solutions were obtained for Reynolds numbers of 100 and 1000. As shown 
in Fig. 8.14, the grids are concentrated near the corners in order to 
capture the secondary vorticities. The results are shown in Figs. 8.15 
and 8.16 for Re^lOO and 1000, respectively. The secondary vorticities 
are captured, but the results are not in good agreement with those 
obtained previously. The maximum value of the primary streamlines is 
30 percent lower than the results obtained with uniform grid distribu- 
tion. Accuracy can be gained by using finer grids but this will result 
in higher computer costs. 












Chapter 9 


CONCLUSIONS 


The main objective of this study has been to investigate the 
feasibility of the method of Lines for application to physical problems 
with nonuniform grid distributions. To attain this objective, it has 
been necessary to investigate the stiffness characteristics of the 
pertinent equations. 

The following conclusions are drawn from the analysis of stiff 
differential equations: (1) equations become stiffer as f^rids are 

concentrated, (2) the slope of the stiffness characteristic curve is 
higher near the fixed boundary (x“0) than near the derivative boundary 
(x»l), (3) equations become stiffer for a large number of grid points, 
(4) use of forward differencing is not feasible, (5) the step size is 
inversely proportional to the Reynolds number for central differenc- 
ing, and (6) backward differencing is preferable over central differ- 
encing at high Reynolds numbers. These conclusions are based on the 
analysis of one-dimensional heat conduction and fluid flow equations. 

The viability and validity of the method of lines are illustrated 
with applications to the Navier-Stokes equations. A computer program 
is developed to solve Eqs. (2.21) and (2.22), and the details of this 
program are available in [ 38] . This method is quite convenient 
from a programming point of view, but computational costs are re- 
latively higher compare to the standard finite difference technique. 

The computational procedure used in this study does not require any 
local linearization in the governing equations. The method has been 
applied successfully to obtain solutions for flow in parallel ducts 


curved-wall diffusers and a driven cavity. Most solutions have been 
obtained by using the non-’conservative form of the governing equations . 

The results are in good agreement with other available results for flow 
in parallel ducts. The solutions also predict the flow development and 
separation in the curved-wall diffusers. For the case of a driven 
cavity, the results obtained by using the non-conser- vative form of 
equations are valid only up to a Reynolds number of 100. At higher 
Reynolds numbers the non-conservative equations provide spurious 
results. The conservative form of the equations was used to obtain 
results up to Reynolds number of 10,000. Both counter rotating 
secondary vorticies are captured by using a nonunlforra grid distribu- 
tion where grids are concentrated near the walls. However, there is a 
loss of accuracy in the primary vortex. The results are in good 
agreement with results of other investigations. 

The analysis of flow in a driven cavity indicates that at a 
Reynolds number of 100, the flow is viscous in the primary vortex with 
little indication of an inviscid core. The inviscid core region has 
filled much of the cavity with viscous regions confined near the walls 
at a Reynolds number of 1000. At a Reynolds number of 10,000, however, 
the inviscid core has completely filled the cavity with small viscous 
regions very close to the wall. These results indicate that the flow 
in a cavity can be characterized by a large vortex adjacent to the 
moving wall and small counter-rotating vorticies in the corners. 

The results show the validity and viability of the method of 
lines where the physical domain is covered with a variable mesh. For 
further study, it is recommended to include the following: (1) use a 
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larger moaber of grid points, (2) solve Che stream function equation by 
a more powerful method, (3) use a higher order approximation for Che 
vorticiCy equation, and (4) write the available codes for a vector 
computer such as Che CYBER 203 so Chat finer grids can be used. 


» 
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